Association of Polygenic Score and the involvement of Cholinergic and Glutamatergic Pathways with Lithium Treatment Response in Patients with Bipolar Disorder

Lithium is regarded as the first-line treatment for bipolar disorder (BD), a severe and disabling mental disorder that affects about 1% of the population worldwide. Nevertheless, lithium is not consistently effective, with only 30% of patients showing a favorable response to treatment. To provide personalized treatment options for bipolar patients, it is essential to identify prediction biomarkers such as polygenic scores. In this study, we developed a polygenic score for lithium treatment response (Li+PGS) in patients with BD. To gain further insights into lithium’s possible molecular mechanism of action, we performed a genome-wide gene-based analysis. Using polygenic score modeling, via methods incorporating Bayesian regression and continuous shrinkage priors, Li+PGS was developed in the International Consortium of Lithium Genetics cohort (ConLi+Gen: N=2,367) and replicated in the combined PsyCourse (N=89) and BipoLife (N=102) studies. The associations of Li+PGS and lithium treatment response — defined in a continuous ALDA scale and a categorical outcome (good response vs. poor response) were tested using regression models, each adjusted for the covariates: age, sex, and the rst four genetic principal components. Statistical significance was determined at P< 0.05. Li+PGS was positively associated with lithium treatment response in the ConLi+Gen cohort, in both the categorical (P=9.8×10–12, R2=1.9%) and continuous (P=6.4×10–9, R2=2.6%) outcomes. Compared to bipolar patients in the 1st decile of the risk distribution, individuals in the 10th decile had 3.47-fold (95%CI: 2.22–5.47) higher odds of responding favorably to lithium. The results were replicated in the independent cohorts for the categorical treatment outcome (P=3.9×10–4, R2=0.9%), but not for the continuous outcome (P=0.13). Gene-based analyses revealed 36 candidate genes that are enriched in biological pathways controlled by glutamate and acetylcholine. Li+PGS may be useful in the development of pharmacogenomic testing strategies by enabling a classification of bipolar patients according to their response to treatment. Keywords: Polygenic score, pharmacogenomics, lithium, bipolar disorder, psychiatry


Introduction
Bipolar disorder (BD) is a severe and often disabling mental health disorder that affects more than 1% of the population worldwide and is characterized by recurrent episodes of depression and mania 1 . BD accounted for 9.3 million disability-adjusted life years (DALYs) in 2017, and imposes a signi cant social and economic burden on society and healthcare systems 2,3 . BD is associated with a signi cant somatic and psychiatric comorbidity 1 and an increased risk of suicide 4 .
Since the discovery of lithium's mood-stabilizing property in 1949 5 , it has been widely used as a rst-line therapy for patients with BD 6, 7 . Lithium is effective in treating acute episodes of illness and reduces the risk of future recurrences of mania and depression 8 . It has also been shown to reduce the risk of suicide 9 . Despite these merits, the e cacy of lithium is highly variable, with about 30% of treated patients showing a favorable response while more than 30% of them have no clinical response at all 8, 10 . Thus far, the causes and predictors of such heterogeneity in treatment response are insu ciently understood.
Genetic factors are thought to contribute, at least in part, to the large interindividual differences in response to lithium [10][11][12][13][14][15] . So far, only a few genetic studies have identi ed speci c single nucleotide polymorphisms (SNPs) and candidate genes associated with patients' response to lithium or treatmentrelated side effects 10,11,[13][14][15][16] . Each employing a genome-wide association study (GWAS) approach, the Taiwan Bipolar Consortium found SNPs in the introns of GADL1 associated with lithium treatment response 17 , whereas the International Consortium on Lithium Genetics (ConLi + Gen) identi ed a locus on chromosome 21 10 , and a follow-up analysis uncovered additional variants within the human leukocyte antigen (HLA) region 14,16 . Gene expression analysis of ConLi + Gen data also showed overexpression of genes involved in mitochondrial functioning in lithium responder patients, highlighting the electron transport chain as a potential target of lithium 18 .
In our recent work, we applied a polygenic score (PGS) modeling approach and demonstrated associations between a poor response to lithium and a high genetic loading for schizophrenia (SCZ) 14 , major depression (MD) 13 , and a meta-PGS combining both SCZ and MD 15 . Machine-learning models that combined clinical variables with the PGS of SCZ and MD has further improved the prediction of lithium treatment response, explaining 13.7% of the variance 19 .
Based on these previous results, translation of PGS testing into clinical practice requires the consideration of three important learnings. First, the PGS of a single phenotype (e.g., SCZ or MD) explains only a small proportion (< 2%) of the variability to treatment response in patients with BD 13, 14 , providing insu cient power for clinical use. Second, a meta-PGS from multiple related phenotypes has better predictive power than a PGS from a single phenotype 15 , suggesting the need to explore additional biological markers, including additional PGSs, that can either independently or together with existing PGSs better predict lithium treatment response. Third, developing polygenic markers with direct pharmacogenomic implications is essential, for example, a PGS for lithium treatment response (Li + PGS ), which is perhaps biologically more related to lithium's pharmacological actions than PGSs built for other clinical phenotypes (i.e. SCZ or MD; that may indirectly in uence treatment response or symptom severity, but do not index pharmacogenetic signatures per se).
Here, we developed a novel Li + PGS for lithium treatment response and applied gene-based pathway analyses to identify molecular mechanisms impacted by genetic variation in response phenotypes.
Findings may assist in optimizing and personalizing the selection of mood stabilizers in patients with BD, and may point to novel molecular targets for future drug development.

Methods And Materials
Study Samples: For this study, we obtained genetic and clinical data from the International Consortium on Lithium Genetics (ConLi + Gen: N = 2,367), Pathomechanisms and Signature in the Longitudinal Course of Psychosis study (PsyCourse: N = 89), and BipoLife cohort (N = 102). Figure 1 shows the detailed steps of data analysis. Insert Fig. 1 Discovery cohort ConLi + Gen is a global collaboration of scientists established to study the pharmacogenomics of lithium treatment in patients with BD 10 . In the current study, we analyzed the genome-wide genotype and clinical data of 2,367 lithium-treated bipolar patients of European ancestry collected by 22 participating sites in 13 countries, including Australia (n = 122), Austria (n = 43), Czech Republic (n = 45), France (n = 210), Germany (n = 218), Italy (n = 255), Poland (n = 97), Romania (n = 152), Spain (n = 74), Sweden (n = 304), Switzerland (n = 57), Canada (n = 353) and the USA (n = 437) 10,20 .

Replication cohort
To replicate Li + PGS associations found in the discovery ConLi + Gen sample, we utilized datasets from PsyCourse and BipoLife where the study participants were of European ancestry. PsyCourse is a longitudinal multicenter study conducted from 2012 to 2019 in Germany and Austria, with up to four assessments at 6 monthly intervals. The study comprises 1,320 patients from psychotic-to-affective spectrum, of which, datasets from 89 patients with BD who received lithium treatment were obtained for this study 21 . BipoLife is a multicenter cohort study, established to investigate the biological basis of BD and patients' response to treatment and being conducted across ten university hospitals in Germany (Berlin, Bochum, Dresden, Frankfurt, Göttingen, Hamburg, Heidelberg, Marburg, Munich and Tübingen) and the medical informatics section of the University of Göttingen 22 .

Target outcome
In both discovery and replication cohorts, patient's treatment response was assessed using the "Retrospective Criteria of Long-Term Treatment Response in Research Subjects with Bipolar Disorder" scale, also called the ALDA scale 10 . The target outcome "lithium treatment response" was de ned in categorical and continuous scales among patients who had received lithium for a minimum of 6 months 10 . The detailed procedures of ALDA scale measurement and its validity are described elsewhere 13,14,20 .
Brie y, the ALDA scale measures symptom improvement over the course of treatment (A-score, range 0 − 10), which is then weighted against ve criteria (B-score that assesses confounding factors, each scored 0, 1, or 2). Once we calculated the total score as 'A-score minus B-score and setting negative scores to zero', the categorical (good versus poor) lithium treatment response was de ned at a cut-off score of 7, where patients with a total score of 7 or higher were considered as "responders" 10 . The continuous outcome for lithium treatment response was de ned on subscale-A, but patients with a total B score greater than 4 or who had missing data on the totals of ALDA subscale-A or B were excluded 10 .

Genotyping, Quality Control And Imputation
We obtained the genotype data assayed with different types of commercial SNP arrays across multiple cohorts 10,21,22 and applied a series of quality control (QC) procedures before and after imputation using PLINK 23 . First, SNPs that had a poor genotyping rate (< 95%), strand ambiguity (A/T and C/G SNPs), a minor allele frequency (MAF) less than 1% or showed deviation from Hardy-Weinberg Equilibrium (P < 10 − 6 ) were removed. Then, individuals with low genotype rates (< 95%), who had sex inconsistencies (between the documented and genotype-derived sex), and who were genetically related were excluded. Imputation: The genotype data passing QC were imputed on the Michigan server 2424 (https://imputationserver.sph.umich.edu) separately for each genotyping platform, using the Haplotype Reference Consortium (HRC) reference panel that consists of the largest available set (64,976 human haplotypes) of broadly European haplotypes at 39,235,157 SNPs 25 . For each cohort, imputation quality procedures were implemented to exclude SNPs of low-frequency (MAF < 10%) and low-quality (imputation quality score R-square < 0.6). From the imputed dosage score, genotype calls for the ltered SNPs were derived and common sets of 4,652,947 SNPs across the cohorts were merged using PLINK 23 .

STATISTICAL ANALYSIS
We implemented polygenic score modeling, genome-wide SNP association, gene-based and functional analysis as described below.
Genome-wide SNP association analysis: Genome-wide SNP association analyses were performed on the binary lithium treatment response and continuous ALDA total score using logistic and linear regression models as implemented in PLINK software 23 , respectively. Each analysis was adjusted for the covariates: age, sex, chip type and the rst four genetic principal components (PCs). For the ConLi + Gen study, Li + PGS was derived only for the European ancestry individuals (n = 2,367) using a ve-fold leave-one-group out (LOG) procedure 28 to remove discovery-target circularity. In each fold, 80% of the sample (n = 1,893) was used to generate GWAS summary statistics that were used as discovery for PGS calculation in the 20% left-out target sample (n = 474). The procedure was repeated ve times by selecting a non-overlapping set of 20% left-out samples to calculate PGS for the entire cohort. Finally, Li + PGS was computed for the PsyCourse and BipoLife participants using ConLi + Gen's GWAS summary statistics (discovery sample) generated from the full European cohort (n = 2,367).
Polygenic score association analysis: To assess the association of Li + PGS with lithium treatment response, a binary logistic regression model was applied for the binary outcome (good versus poor response to lithium treatment), and a Tobit analysis model (censored regression) was used for the continuous outcome (ALDA total) 29 . In addition, we divided the ConLiGen sample into deciles, ranging from the lowest polygenic load (1st decile, reference group) to the highest polygenic load (10th decile).
Then, we compared BD patients in the higher polygenic load deciles (2nd − 10th deciles) with patients in the lowest polygenic load decile (1st decile). In both the binary and continuous outcomes, the proportion of phenotypic variance explained by Li + PGS was computed as the difference in R 2 of the model t with Li + PGS plus covariates, compared to the model t with only covariates. Each modeling analysis was adjusted for the covariates: age, sex, and the rst four genetic PCs, and statistical signi cance was set at p < 0.05.

Gene-based and functional analysis
The gene-based analysis was based on summary statistics generated from the full European ancestry (n = 2,367) ConLi + Gen genome-wide SNP association analysis (see Fig. 3) and employed MAGMA (Multimarker Analysis of GenoMic Annotation) 30 , a tool that uses a multiple regression approach to incorporate LD between markers and to detect multi-marker effects.
To explore the biological context of the genes discovered from the gene-based analysis, a pathway analysis was implemented using PANTHER (Protein ANalysis THrough Evolutionary Relationships; http://pantherdb.org/) classi cation system. PANTHER is designed to classify proteins (and their genes) into biological pathways 31 . To prepare the input genes for PANTHER, we selected genes that showed gene-level association with lithium treatment response (either with the categorical or continuous outcome) at MAGMA adjusted p-value < 0.001. This list of genes was entered into PANTHER version-17 which compares the proportion of input genes mapping to a biological pathway to the reference gene list from its databases. Molecular relationships previously experimentally observed in Homo sapiens (human) were included. The signi cance of the overrepresented PANTHER pathways was determined using Fisher's exact test and later adjusted for multiple testing using the Bonferroni correction method.
Signi cant associations were de ned at p-value < 0.05.

Sample Characteristics
The discovery analysis consisted of ConLi + Gen data obtained from 2,367 bipolar patients of European ancestry who had undergone lithium treatment for at least six months. The mean (sd) age of the patients was 47.5(13.9) years and 1,369 (57.8%) were female. In all, 660 (27.9%) of patients had a good response Page 14/29 to lithium treatment (ALDA score ≥ 7). The mean (sd) ALDA score for ConLi + Gen participants was 4.1 (3.1). The replication analysis was based on a combination of the PsyCourse and BipoLife datasets (N = 191), whose mean (sd) age was 49.1(13.0) years. Of the 191 patients with BD, 48(25.1%) had a good response to lithium (Table 1).  (Table 2). While there was an increasing trend in the odds of lithium treatment response across the deciles, the most signi cant prediction contrast was found at the 'extremes' (1st and 10th decile) which comprised of ~ 20% of the total cohort (Fig. 2). A replication PGS analysis in the combined PsyCourse and BipoLife samples found a statistically signi cant association of Li + PGS with the categorical lithium treatment response (P = 3.9x10 − 4 , R 2 = 0.9%), but not with the continuous outcome (P = 0.13).  Fig. 2 here Genome-wide association, gene-based and functional analysis After re-imputing the ConLi + Gen data in reference to the latest HRC genomes, we conducted GWASs on lithium response, both in categorical and continuous outcomes. This GWAS analysis identi ed a single locus with lead SNP rs9396756 located near the stathmin domain containing 1 (STMND1) gene that reached genome-wide signi cance for association with the categorical outcome (P = 2.7 x10 − 8 ) and showed a suggestive association with the continuous ALDA score (P = 7.6 x10 − 8 ) (Fig. 3). A follow-up gene-based analysis of the newly derived ConLi + Gen GWAS summary statistics found 36 candidate genes likely associated with lithium treatment response -assessed in either continuous or categorical outcomes (P < 0.001). In silico functional analysis of the 36 genes revealed enriched biological pathways including the muscarinic acetylcholine receptors 1 and 3 (P-value corrected for multiple testing = 0.026) and metabotropic glutamate receptor group III pathway (P = 0.043). These genes and pathways may have an impact on clinical response to lithium treatment and be potential molecular targets for lithium (Supplementary Fig. 1

Discussion
This study presents ndings from a comprehensive analysis of genetic and clinical data on lithium treatment response that involved the development of a polygenic score for lithium treatment response (Li + PGS ), genome-wide SNP association and gene-based and functional analyses.
Since the publication of the rst GWAS report by the ConLi + Gen team 10 , two landmark studies that independently showed the negative association of PGSs for SCZ and MD with lithium treatment response have been published [13][14][15] . The rst study found that 10% of bipolar patients with the lowest polygenic load for SCZ were 3.46 times more responsive to lithium compared to 10% of patients with the highest genetic load for SCZ 14,15 . Similarly, in the second study, 10% of patients who had the lowest genetic loading for MD were 1.54 times more responsive to lithium than 10% of patients with the highest genetic loading for MD 13,15 . Nevertheless, each of these PGSs accounts for < 2% of the total variance to lithium treatment response 13 , suggesting the need to explore additional biological traits that can either independently, or in concert with existing PGSs better predict lithium response. Moreover, the previous PGSs for SCZ and MD are di cult to interpret in a pharmacogenomic context, making the development of a speci c lithium response PGS necessary, which is assumed to be more likely to be associated with lithium treatment response and perhaps is biologically more related to lithium's pharmacological actions.
In this novel study, we constructed a PGS for lithium response-Li + PGS, a biological marker of direct pharmacogenomic relevance, and showed a positive relationship between a high genetic loading for lithium treatment response variants and long-term therapeutic response to lithium in patients with BD. We demonstrated that bipolar patients at the extreme tail end of the distribution have the strongest association, i.e. 10% of patients who carry high genetic loading for lithium responsive variants (10th decile) were 3.47 times more likely to respond to lithium compared to 10% of those with the lowest Li + PGS (1st decile). These results indicated that Li + PGS has the potential to help stratify bipolar patients according to predicted lithium response.
In a GWAS of lithium treatment response, we identi ed a locus near the STMND1 gene, which encodes for proteins known to be involved in neuron projection development, and active in neuron junctions and cytoplasm. Previous analysis that employed the 1000 Genomes Project reference panel for imputation reported a suggestive association between genetic variants within the STMND1 gene and lithium treatment response 10 .
Using our newly derived ConLi + Gen GWASs summary statistics as an input, we then carried out a genebased analysis where several genetic variations were examined together for their association with lithium treatment response 30 . This approach found 36 potential target genes for lithium treatment that are enriched in the muscarinic acetylcholine receptors (mAChRs) 1 and 3 and the metabotropic glutamate receptor group III signaling pathways -well characterized biological pathways modulated by the most abundant neurotransmitters in the brain (glutamate and acetylcholine).
Acetylcholine is the central regulator of the mAChRs signaling pathways, which are subfamily of G protein-coupled receptor complexes located in the cell membranes of neurons and other cells that regulate fundamental functions of the central and peripheral nervous system including acting as the main end-receptor stimulated by acetylcholine released from postganglionic bers in the parasympathetic nervous system 32 . The muscarinic antagonist scopolamine has antidepressant activity, while physostigmine, a cholinesterase inhibitor induces depressive symptoms, suggesting muscarinic receptors may play a role, not only in the pathogenesis of mood disorders, but also as therapeutic targets 33 . M1 and 3 receptors are localized in the cortex, hippocampus and substantia nigra and are known to activate protein kinase C (PKC), causing post-synaptic excitation. PKC is thought to be central in the molecular pathogenesis of BD.
Glutamate, the primary excitatory neurotransmitter in the central nervous system (CNS), exerts neuromodulatory actions via the activation of metabotropic glutamate (mGlu), a type of glutamate receptor that modulates synaptic transmission and neuronal excitability throughout the central nervous system 34 . Group III metabotropic glutamate receptors are largely presynaptically localized and downregulate neurotransmitter release from presynaptic terminals directly or indirectly. These receptors have a prominent expression in the brain, especially in the region of the hippocampus, and can lead to the inhibition of the cAMP cascade which is critical for the maintenance of long-term synaptic plasticity 35 .
Growing evidence indicates that abnormalities in the glutamatergic system are implicated in the pathogenesis and treatment of mental health disorders 36 including BD 37,38 , SCZ 39 , neurodevelopmental disorders 40 , Huntington's disease 41 and Alzheimer's disease 42 . Studies have reported SNPs of the mGluRs system associated with BD 43 , and in animal studies, lithium was found to alter intracellular calcium by modulating the activity of the metabotropic glutamatergic receptor system 44 . To summarise, ndings from the genome-wide SNP association, gene-based and functional analysis highlight the possibility that mechanisms involving glutamate and acetylcholine signaling pathways might in uence the therapeutic effects of lithium in patients with BD. Modulation of these pathways through genetic variants may disrupt or enhance lithium's clinical effectiveness.
Our study has some limitations. First, while our ndings were replicated in an independent small size sample, the fact that it was replicated in the binary outcome, but not in the continuous outcome indicates the need for a larger replication cohort. Second, because Li + PGS was developed and evaluated in also received a research grant from Takeda Pharmaceutical Co, Ltd. Peter Falkai has received grants and served as consultant, advisor or CME speaker for the following entities Abbott, GlaxoSmithKline, Janssen, Essex, Lundbeck, Otsuka, Gedeon Richter, Servier and Takeda as well as the German Ministry of Science and the German Ministry of Health. Eva Reininghaus has received grants and served as consultant, advisor or CME speaker for the following entities: Janssen and Institut Allergosan. Mikael Landén declares that, over the past 36 months, he has received lecture honoraria from Lundbeck and served as a scienti c consultant for EPID Research Oy; no other equity ownership, pro t-sharing agreements, royalties or patent. Kazufumi Akiyama has received consulting honoraria from Taisho Toyama Pharmaceutical Co, Ltd. In 2021, Jörg Zimmermann served as an advisor for Biogen concerning Aducanumab (Alzheimer`s Disease).The other authors have no other con ict of interest to disclose.
LOG=Leave-one-group out procedure; PsyCourse = Pathomechanisms and Signature in the Longitudinal Course of Psychosis study and BipoLife = German research consortium for the study of bipolar disorder.

Figure 2
Trends in the odds ratios (ORs) for favourable treatment response to lithium for patients with bipolar disorder in the higher genetic loading for lithium responsive variants, deciles (2 nd to 10 th ) compared with patients in the lowest (decile 1 st ) of genetic loading for lithium response (n = 2,367).
Legend: The X mark on the line plot indicates that the association is not statistically signi cant at that decile.